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ABSTRACT 

We apply a semianalytic galaxy formation model to two high resolution cosmological N-body 
simulations to investigate analogues of the Milky Way system. We select these according to 
observed properties of the Milky Way rather than by halo mass as in most previous work. 
For disk-dominated central galaxies with stellar mass (5-7) x 1O 1O M 0 , the median host halo 
mass is 1.4 x 1O 12 M 0 , with la dispersion in the range [0.86, 3.1] x 1O 12 M 0 , consistent with 
dynamical measurements of the Milky Way halo mass. For any given halo mass, the proba¬ 
bility of hosting a Milky Way system is low, with a maximum of ~20% in haloes of mass 
~ 1O 12 M 0 . The model reproduces the F-band luminosity function and radial profile of the 
bright (My < —9) Milky Way satellites. Galaxy formation in low mass haloes is found to be 
highly stochastic, resulting in an extremely large scatter in the relation between My (or stellar 
mass) for satellites and the depth of the subhalo potential well in which they live, as measured 
by the maximum of the rotation curve, F max . We conclude that the “too big to fail” problem 
is an artifact of selecting satellites in iV-body simulations according to subhalo properties: 
in 10% of cases we find that three or fewer of the brightest (or most massive) satellites have 
F max > 30kms _1 . Our model predicts that around half of the dark matter subhaloes with 
F max > 20 km s -1 host satellites fainter than My = — 9 and so may be missing from existing 
surveys. 

Key words: 


1 INTRODUCTION 


The cold dark matter (CDM) model has been shown to be consis¬ 
tent with many observations on cosmological scales, but uncertain¬ 
ties remain on small scales where complex baryonic processes are 
involved. A substantial number of faint satellite galaxies are known 
in the Milky Way system. Thanks to their relative proximity, the 
distribution, chemistry and motions of their individual stars can be 
measured precisely. This makes the Milky Way satellites an excel¬ 
lent laboratory to test the CDM model and also to investigate the 
physics of galaxy formation on small scales. 


It has been known for some time that the count of substruc¬ 
tures in simulated Milky Way-mass dark matter haloes greatly ex¬ 
ceeds that of the known luminous Milky Way satellites (e.g. |Klypin| 
|et al. 1 1999) [Moore et al.|1999} . In the context of CDM this is read¬ 
ily explained if the efficiency with which baryons are converted 
into stars drops quickly at low halo masses. This inefficiency is 
expected on the basis of the well understood atomic physics gov¬ 
erning radiative cooling, the existence of an ionising cosmic UV 
background, and comparison of the energy released by supernovae 
(which drive the expulsion of baryons from dark matter haloes) to 
the depth of halo potential wells i Efstathiouj 1992, KauffmannetalJ 


[19931 |Bullock et al.|2000| |Benson et al.|2002| |Somerville|2002| >. 


CDM galaxy formation models which incorporate even these most 
basic astrophysical effects have been able to reproduce not only the 
abundance but also radial distribution of satellites around the Milky 


Way (e.g. 

Okamoto et al.J2010, Maccio et al. 2010] [Li et al. 2010; 

Guo et al. 

2011, Font et al. 2011, Parry et al. 2012). A host of other 


effects, such as cosmic ray pressure, may contribute to the regula¬ 


tion of star formation in these most marginal systems, (|Wadepuhl| 
|& Springel|2011) . 


Dwarf spheroidal galaxies (dSphs) make up most of the 
Milky Way (MW) satellite population. Since they are dark matter- 
dominated galaxies, they are ideal for testing the close connection 
between the properties of dark matter haloes and the assembly of 
stellar mass predicted by CDM models. Relatively direct compari¬ 
son between models and observations is now possible thanks to de¬ 
tailed kinematic analyses of the Milky Way dSphs (e.g. Penarrubia 


et al.|2008||Strigari et al.|[2008l |Lokas|2009l |Walker et al.|2009 


Wolf et al.|2010||Strigari et al.|2010| . |Boylan-Kolchin et al. 1(201 I 
201 2[ > compared observed dSph stellar kinematics to predictions 
from the Aquarius Project, a set of six N-body simulations of dark 
matter haloes of mass ~ 10 12 Mg (consistent with constraints on 
the halo mass of the Milky Way; see e.g. figure 1 of Wang et al. 
2014). They concluded that the most massive subhaloes in such 
simulations, which in typical galaxy formation models would host 
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the most luminous satellites, are always too dense to be dynami¬ 
cally consistent with observations of any of the known Milky Way 
companions. [Boylan-Kolchin et al.|(2011[ l dubbed this discrepancy 
the “too big to fail” (TBTF) problem. 

Several possible solutions to the problem have been advanced, 


including alternative forms of dark matter (e.g.|Lovell et al. |2012 | 
Vogelsberger & Zavala|20 1 3), baryon-induced changes in dark halo 


density profiles (jdi Cintio et al. 2011 

Garrison-Kimmel et al. 2013 

Arraki et al. [2013, Brooks & Zolotov 

20141, and uncertainties in the 

mass of the Milky Way halo (Wang et al. 2012 [ Cautun et al. 2014). 


All of these studies use dark matter halo mass as a starting point to 
select Milky Way analogues in simulations in order to analyze the 
internal kinematics of their subhaloes. However, galaxy formation 
involves many complex processes and its outcome depends not only 
on the present-day halo mass, but also on the formation history of 
the system. This is an important consideration both for the selection 
of primary galaxies and for the identification of relevant satellite 
haloes. Indeed, [Sawala et al.| < j2014) t show that huge scatter exists 
in the relation between the stellar mass of dwarf satellite galaxies 
and their present-day subhalo mass. |Sawala et aL] ( |2014| ) simulated 
analogues of the Local Group including full baryonic physics and 
found that the expulsion of baryons in winds from small dark matter 
haloes at early times reduced the central density of the haloes in 
which satellites form enough to solve the TBTF fail problem. 

In this paper we implement the galaxy formation model of 
|Guo et ah] < )20 13) > on two N-body cosmological simulations and 
identify Milky Way analogues using the properties of their cen¬ 
tral galaxies. We compare the halo mass of these Milky Way ana¬ 
logues to the observations, and examine the possibility of dark mat¬ 
ter haloes of given mass hosting the Milky Way analogues. We then 
further analyze the abundance and kinematics of the subhaloes of 
satellite galaxies selected by their luminosity or stellar mass. At the 
end we investigate whether selecting MW analogues according to 
these observables could help solve the TBTF problem. 


2 SIMULATION AND SEMI-ANALYTIC MODELS 

This work makes use of two ACDM simulations: a high-resolution 
zoom simulation, Copernicus Complexio (COCO; Hellwing et al. 
in prep), and a cosmological volume simulation of lower resolu¬ 
tion, Millennium-II (MS-II [Boylan-Kolchin et al.[2009) . The MS- 
II cube has sidelength 100/i _i Mpc. The particle mass is m p = 
6.9 x 10 6 /i _1 Mq and the force softening scale is e = l/i -1 kpc. 
It assumes WMAP-1 cosmological parameters, with a linear power 
spectrum normalization, <T8 = 0.9, and matter density, = 0.25. 

COCO simulates a high resolution comoving volume V = 
2.25 x 10 4 /i _s Mpc 4 in the centre of a lower resolution peri¬ 
odic cosmological cube of sidelength 70.4/i _1 Mpc using N p = 
2374 s DM particles. This results in a sample of ~ 50 haloes of 
~ 10 12 Mg at a = 0 with high mass (m p = 1.135 x 10 5 1 i _1 Mq ) 
and force (e = 230/U 1 pc) resolution. Initial density perturbations 
were generated with the novel Panphasia technique ( |Jenkins|20 1 
which provides self-consistent and reproducible random phases for 
zoom-in resimulations across an arbitrary range of resolutiorQ 

Unlike MS-II, COCO assumes WMAP-7 cosmological pa¬ 
rameters (O m = 0.272, Oa = 0.728, as = 0.81 andn s = 0.968). 
The slightly different cosmogonies simulated by these two calcu¬ 
lations are reflected in the abundances and internal properties of 

1 The specific COCO phase information is available on a request from the 
authors. 


their dark matter haloes. For example, the abundance of MW-mass 
haloes differs by a few per cent ( jBoylan-Kolchin et al.|201 1~[ . How¬ 
ever, the structural and kinematic properties of massive subhaloes 
in MW haloes are not affected in any significant way jB oylan- 

Kolchin^taH20Ll| |Wang et al.|2012||Garrison-Kimmel et al.|20T3 

Cautun et al.|2014). 


We used stored snapshots (160 for COCO and 68 for MS- 
II) for Friend-of-Friends group finding jDavis et al. 1 1985) with a 
linking length equal to 0.2 times the mean inter-particle separa¬ 
tion. Subhaloes are identified by applying the SUBFIND (Springel 
|et al.pOOl) algorithm to each FOF group. We define the centre 
of the FoF group as the potential minimum of the most massive 
self bound sub halo. The viral mass, M v i r , is defined as the mass 
enclosed by a virial radius, R v i r , which we approximate by the ra¬ 
dius, R 200 , within which the mean density is 200 times the crit¬ 
ical value for closure. The maximum circular velocity of a halo 
is defined as Unax = max [^/GM (r) /r] and is attained at radius 
v max - For haloes close to the resolution limit, V ma , x can be systemat¬ 
ically underestimated. To correct for this effect, following|Springel| 


et al. (2008), we adjust measurements Umax from the simulation to 

Umax = u^ ax (l + (e/r max ) 2 ) 0 ' 5 . 


To populate dark matter haloes in our two simulations with 
galaxies, we applied the semi-analytic model recently developed 
by the Munich Group JGuo et aT[|20 1 1 1 12013) to their merger 
trees. This model reproduces many statistical properties of the ob¬ 
served galaxy population, including the abundance of bright satel¬ 
lites around the Milky Way l |Guo et al.|201 1) . There are two pro¬ 
cesses crucial in understanding the formation of low mass systems: 
UV reionization and supernova (SN) feedback. |Guo et al. |201l| 
adopted a fitting function originally proposed by |Gnedin| 


2013 


et al. 


(|2004[> to describe the baryon fraction as a function of halo 


mass and redshift. The fitting parameter and the characteristic halo 
mass beyond which the baryon fraction is close the universal value 
but below which it rapidly drops with decreasing halo mass are 
given by |Okamoto et al.|(2008} . 


|Guo et ah] ( |2011| |2013| > introduced a SN feedback model 
which depends strongly on the maximum velocity of the subhalo, 
but the total amount of energy to reheat and eject gas can never 
exceed the total energy released by SNII. For most of the satel¬ 
lites around the Milky Way, the feedback energy saturates at this 
value. Guo et al.] ([201 lj also showed that SN feedback dominates 
the formation of relative luminous satellites (Mv < -11), while 
reionization becomes important only for the very faint satellites 
(Mv > -ID. 


In |Guo et al.| ( |2011) , satellite galaxies within the virial radius 
of their host are subject to two environmental effects: stripping of 
their hot gas halo by tides and ram pressure, and tidal disruption of 
their stellar component. The orbit of the satellite’s subhalo in the 
A-body simulation is used to estimate these tidal and ram pressure 
forces. Once a satellite galaxy loses its parent subhalo, it is consid¬ 
ered to merge into the central galaxy unless the density of the host 
halo at the pericenter of the satellite’s orbit exceeds the average 
baryonic mass density within the half mass radius of the satellite. 
In this latter case, cold gas from the disrupted satellite is added to 
the hot gas halo of the host, and stars from the satellite to a ‘stellar 
halo' component (which is not counted towards the stellar mass of 
the central galaxy). 


In the following sections we describe how we select MW ana¬ 
logues from the resulting mock catalogues and present our analysis 
of their satellite galaxy populations. 
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3 RESULTS 

The particle mass resolution of MS-II is sufficient to study the for¬ 
mation history of haloes more massive than ~ 10 10 Mq. More¬ 
over, in our semi-analytic model, the properties of central galax¬ 
ies in haloes of mass 10 12 Mg are converged at the resolution 
of MS-II. However, the resolution of MS-II is not sufficient to 
study the internal kinematics of satellites with Umax < 30 kins -1 
( |Wang et al.|2012} . The resolution limit of COCO corresponds to 
Knax ~ 10 kms -1 (Hellwing et al. in prep), corresponding to 
the least massive satellite galaxies known in the Milky Way sys¬ 
tem. The smaller volume of COCO yields relatively few Milky Way 
analogues, however. We therefore use the larger volume of MS-II to 
study the halo mass distribution of Milky Way analogues selected 
according to their observable properties, and the higher resolution 
of COCO to study the internal properties of subhaloes in systems 
selected in this way. 


3.1 Halo mass of the Milky Way 


In previous works Milky Way analogues have been identified by 
selecting isolated haloes with mass in a narrow range, typic ally 
around ~ 10 12 Mq. However, both theoretical work (e.g. Guo 


et al.|2010l|Moster et al.|2010||Behroozi et al. 2010 Cautu n et al. 


2014[> and observa tional data (e.g. Mandelb aum et a ,|2006|[Leau 
thaud et al.|2012| > suggest that the scatter in the relation between 
galaxy mass and host halo mass is large. This raises the question of 
whether or not a narrow halo mass range is sufficiently representa¬ 
tive of the ‘Milky Way’ galaxy population. 

In the following we use galaxy properties predicted by our 
semi-analytic model to select Milky Way analogues, rather than 
dark matter halo mass. Our selection is based the stellar mass and 
morphology of the central galaxy with the following criteria: 


5 x 1O 1O M 0 < M, < 7 x 1O 1O M 0 (1) 

and 

0.05 < M bulge /M» < 0.4, (2) 

where M* is the total stellar mass of the galaxy and M bu i ge 
is the stellar mass of the galactic bulge. We have chosen this stellar 
mass range based on recent observational constraints (e.g. |McMil-| 
lan |201 l| l. The bulge selection is only intended to restrict the sam¬ 
ple to disk dominated galaxies. The actual bulge mass fraction of 
the Milky Way is uncertain, as is the extent to which the model 
‘bulge’ mass corresponds to the photometrically or kinetically clas¬ 
sified ‘bulges’ of real disk-dominated galaxies. This joint selection 
yields 566 Milky Way-like galaxies in MS-II and 13 in COCO. 

In Fig. [T] the middle panel shows the distribution of M v i r for 
simulated Milky Way analogues selected according to these crite¬ 
ria. The black curve represents the MS-II sample and the red curve 
the COCO sample. Although the mass and spatial resolution of the 
two simulations differ greatly, a K-S test supports the hypothesis 
that these two samples are drawn from a common parent distribu¬ 
tion at a confidence level of 0.97 (roughly within Poisson errors, as 
shown). 

The median halo mass is 1.4 x 10 12 Mq, with the 16 th to 
84 th percentile range (8.6x 10 11 Af 0 to 3.1x10 12 Mq) indicated 
by black dashed lines. A number of observational measurements 
of M 200 are shown by symbols with error bars in the top panel of 
Fig. □ These measurements have been made using the kinematics 
of different halo tracers by|Xue et al.|(2008 black circle), [Gnedinj 


et al. (2010 blue upward triangle). Watkins et al. (2010 blue down- 

ward triangle), Sakamoto et al. (2003 

cyan diamonds, with or with- 

out the inclusion of Leo I) and Battag 

ia et al. ( 20051 purple squares. 


assuming a truncated flat model or NFW model); using the escape 
velocity of the Large Magellanic Cloud or satelli tes (|Busha et al.| 
|201 1 1 red diamond); using masers by Klypin et al. (2002 green (' 


cle), incorporating photometric and kinematic data by (McMillan 


2011, green upward triangle); using abundance matching by Guo 


et al.lpOlOl cyan downward triangle); using the escape velocity of 


halo stars by |Piffl et al.| l |2014] black rightfacing triangle); and us¬ 
ing the timing argument of the Milky Way/Andromeda pair by|Li| 
|& White| (2008 red square). Almost all of the measurements in the 
literature fall between the 16 th and 84 th percentiles of the distri¬ 
bution predicted by our model. The model predicts that haloes of 
even higher masses could also host galaxies satisfying our Milky 
Way selection criteria: 31% of the host haloes in our samples are 
more massive than 2 x 10 12 Mq and 17% are more massive than 
3 x 10 12 M q . 

We have tested how a systematic shift in the stellar mass of 
the Milky Way affects the resulting halo mass distribution by as¬ 
suming an alternative stellar mass range, 4 x 10 10 Mq < M* < 
6 x 1O 1O M0. Relative to the peak of the black curve, we find the 
distribution of M vb shifts to lower mass by ~ 0.1 dex. In other 
words, the host halo mass range does not change much if we re¬ 
duce the median ‘Milky Way’ stellar mass by 20%. 

The bottom panel of Fig.^shows the fraction of all dark mat¬ 
ter haloes of a given mass that host a central galaxy satisfying our 
Milky Way selection criteria. We find that this fraction reaches a 
maximum of around 16% at ~ 2 x 10 12 Mq. Most of the haloes 
this mass either host a galaxy substantially more or less massive 
than the Milky Way, or else host a bulge-dominated galaxy. The 
fraction of Milky Ways drops rapidly at lower halo masses, and 
more slowly at higher masses. In general, even in the most favoured 
range of M 200 , the probability for a halo of any given mass to host 
a MW analogue is remarkably low. Hydrodynamic simulations of 
disk galaxies often preselect haloes from cosmological volumes for 
‘zoom’ resimulations studies based on their mass. Our results indi¬ 
cate that such a selection will always be inefficient. Clearly, other 
properties of haloes and their individual formation histories also 
need to be taken into account jSawala et al.|2014). 


3.2 Satellite galaxies around the Milky Way 

In the last section we have shown that our model predicts a typi¬ 
cal halo mass for Milky Way analogues in the same range as re¬ 
cent measurements of the Milky Way galaxy itself. In this sec¬ 
tion, we will focus on the properties and spatial distribution of 
satellite galaxies around the Milky Way. As mentioned in Sec. 2, 
the mass resolution of the COCO simulation is 60 times higher 
than that of MS-II, and the spatial resolution is higher by a fac¬ 
tor of 4. More than 90% of satellites in Mikly Way systems with 
Fmax > 15 km s’ 1 still have their own subhaloes in COCO, which 
allows us to follow their evolution in detail (see Appendix). This 
level of resolution is especially important for studying the internal 
dynamics and structure of satellite galaxies. For example, their ro¬ 
tation curves can be traced reliably and hence I4 lax can be used as 
a robust measure of the depth of their potential. Therefore, in the 
following sections, we only use the COCO simulation to study the 
properties of satellite systems. 
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Figure 1. Upper panel: collection of measurements of M v i r for the Milky 
Way taken from the literature (see text for details). Middle panel: the distri¬ 
bution of host halo mass for model galaxies with 5 < M* < 7 x 10 10 Mq 
and 0.05 < Mb u i ge /Af* < 0.4 selected in the MS-II (black curves) and 
COCO (red curves) simulations. Error bars show the Poisson errors for each 
mass bin. A KS test supports the hypothesis of a common parent distribution 
in the MS-II and COCO samples at a confidence level of 0.97. Dotted lines 
mark the 16 th to 84 th percentile range of the MS-II sample. Bottom panel: 
the probability for a halo of a given mass to host a galaxy that matches our 
Milky Way criteria in MS-II. 


3.2.1 Abundance and profile 

Fig. 0 shows cumulative counts of satellite galaxies as a function 
of V-band magnitude around the 13 simulated MW analogues se¬ 
lected from COCO (black curves). The red dot-dashed curve repre¬ 
sents the equivalent luminosity function of known satellites around 
the Milky Way and the dashed curve shows the average luminosity 
function of satellites in the MW and M31 found by |Koposov et al.| 
< |2008j >, who attempted to correct for incompletness and partial sky 
coverage. For this comparison, following Koposov et al., we de¬ 
fine all galaxies within 280 kpc of each MW analogue as its satel¬ 
lites. The model predictions are broadly consistent with the data 
although the slope at the faint end is steeper in the simulations than 
in the data. This demonstrates convergence with the results of|Guo| 
|et al.|p0TT) , who showed that the bright end of the V-band satllite 
luminosity function for the much larger number of MW analogues 
in MS-II is also in reasonable agreement with the MW and M31 
systems'^] 

FigJ3] compares the galactocentric radial distribution of the 

2 The convergence of the MS-II MW satellite luminosity functions, 
demonstrated here by comparison to COCO, suggests that the semi-analytic 
treatment of orbits for galaxies whose dark matter haloes are stripped below 
the resolution of the IV-body simulation works reasonably well. 



Figure 2. Cumulative distribution of V-band magnitude for satellite galax¬ 
ies in each simulated MW system (black curves). The red dot-dashed line 
shows the equivalent distribution for known Milky Way satellites and the 
red dashed line the average distribution for Milky Way and M31 satellites 
corrected for completeness and sky coverage by |Koposov et al.|f2008f 


most luminous satellites around our Milky Way analogues to those 
in the real Milky Way system. We restrict this comparison to the 
‘classical’ brightest 12 satellites, using the galactocentric distances 
given by |McConnachie| (20 1 2} . These observed satellites galaxies 
are brighter than My ~ —9, so we select simulated satellite galax¬ 
ies with My < —9 for this comparison. The Milky Way obser¬ 
vations lie within the envelope of the radial distributions from our 
model. 


3.2.2 Distribution ofV ma , x 

In the following we focus on the kinematic properties of satel¬ 
lite galaxies. Recently, detailed kinematic measurements have been 
carried out for dSphs in the Local Group (e.g PenamibiaetaL|2008| 


Strigari et alpM) |Lokas||2009| |Walker et al.||2009| |Wolf eTaT] 


2010j |Striganetal 4 20!0| |2014| l, which makes it possible to per¬ 

form a relatively direct comparison between data and simulations. 
In particular, stellar kinematics are often used to infer Knax, which 
is then used to estimate the gravitational potentials of satellite host 
dark matter haloes. 

Vmax is not a direct observable, of course, and must be in¬ 
ferred from other properties, such as the mass within a given radius 
(usually the mass within the half-mass radius, JW 1 / 2 ), which in turn 
is computed from a measured stellar velocity dispersion. To cali¬ 
brate this procedure, Boylan-Kol chin et al.| (2012l computed I4iax 
using subhaloes from the Aquarius suite, a series of very high reso¬ 
lution resimulations of dark matter haloes with M200 ~ 10 12 Mg. 
Assuming that the simulated subhaloes together constitute a repre¬ 
sentative sample of Milky Way-like satellite galaxy hosts, they con¬ 
structed a theoretical distribution of I/ max for each real MW satel¬ 
lite. This was done by assigning a weight to the Knax of each sim- 
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Figure 3. The radial profile of the 12 brightest Milky Way satellite galax¬ 
ies (red dashed curve, McConnachie et al. 2012) compared to the profiles 
of satellites with My < — 9 in our COCO Milky Way analogue systems 
(black curves). 


Figure 4. Cumulative distribution of maximum circular velocity for the 12 
most massive satellite galaxies in each simulated MW system ranked by 
stellar mass (blue and black curves). Red points show the values of Vmax 
for MW satellites given by Boylan-Kochin et al. (2012). The vertical error 
bars on these points show the Poisson error. 


ulated subhalo, according to how closely its mass within a radius 
equal to the observed stellar half-mass radius matched the value of 


Mi/2 inferred from kinematic observations (for details seelBoylan- 
|Kolchin et al.|20l2) . 

Fig. El shows the cumulative distribution of V max for the 12 
classical brightest satellites around the Milky taken from Table 2 
of |Boy lan-Kolchin et al. (2012, red symbols with error bars)). For 
comparison to this data, we rank the most massive 12 satellite 
galaxies of each model Milky Way system according to their stellar 
mass and construct Knax distributions (black and blue lines). The 
lowest value of Knax corresponding to resolved satellites in COCO 
is about 10 kms” (see Appendix). 

Most of the simulated examples overpredict the abundance of 
satellites with high Knax relative to the Milky Way. At K max ~ 10- 
15kms -1 , predicted subhalo abundances are consistent with the 
MW system. At Knax ~ 17 kms -1 , several model examples still 
agree with the data within Poisson errors. One out of the 13 simu¬ 
lated Milky Way analogues has a Knax distribution consistent with 
the data (within Poisson errors) over the whole magnitude range 
considered here. We mark this system with a blue curve in Fig. [4] 

The host halo mass of the satellite system represented by the 
blue curve is 9.7 x 10 11 Mq, around the median of current estimates 
of the MW host halo masses (Fig. [TJl. Flence, this is not the least 
massive halo in our sample. Of all our examples, this system has 
the lowest abundance of satellites with high Knax- This may seem 
somewhat at odds with the conclusions of| Wang et al.| (|2012)l and 
|Cautun et al.| ( {2014] >, who found that the abundance of subhaloes at 
high Knax is proportional to the host halo mass. Although such a 
relation holds on average, several factors introduce a large amount 
of scatter into the correlation. In particular, the luminosity of dwarf 


galaxies in our semi-analytic model depends strongly on their for¬ 
mation history and not just on their present-day halo mass (e.g. 
|Li et al.|2010| l. This broadens the relation between luminosity and 
halo mass and hence that between luminosity and Knax- A similar 
scatter in the properties of low mass galaxies at a fixed halo mass 
is also found in recent hydrodynamical simulations (jSawala et al.J 
|20l4l >. 

3.2.3 Too big to fail? 

|Boylan-Kolchin et al.|j20 1 1) compared the abundance of dark mat¬ 
ter subhaloes as a function of K ma x in the Aquarius simulations 
(assumed to be a representative sample of MW host haloes) to 
the inferred Knax distribution of Milky Way satellites. They found 
that Aquarius predicts significantly more subhaloes with high K m a X 
compared to what one could expect from all-sky extrapolation 
of the available Milky Way data. Specifically, there are only 3 
known satellites galaxies with K max above 30kms -1 (the Large 
and Small Magellanic Clouds, and the tidally disrupted galaxy in 
Sagittarius), whereas Aquarius predicts, on average, 8 subhaloes 
with Knax > 30kms -1 . This is one of a number of ways to ex¬ 
press the “too big to fail” (TBTF) problem mentioned in the intro¬ 
duction. 

|Boylan-Kolchin et al.| < |2012j l restated the TBTF problem by 
comparing the circular velocity profiles of simulated subhaloes to 
the kinematically best-matching observed MW dwarf spheroidal. 
Since the instantaneous value of Knax for a given subhalo evolves 
over time (increasing as a subhalo grows, and decreasing as it is 
tidally stripped), they compare with values Knax(Aref) at several 
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Figure 5. Cumulative distribution of maximum circular velocity for the 
12 satellite galaxies with the highest Knax (upper panel) and stellar mass 
(lower panel) in each simulated MW system, averaged over all simulated 
MW systems. Different colours are for different K ma x cuts as labelled in 
the right bottom comer in the bottom panel. Dotted vertical lines indicate 
the ‘critical number’ of three satellites. 


different reference redshifts (2 re f = 0,10, Zinfaii, etc.). The idea 
is that the present-day stellar mass in a subhalo should be tightly 
correlated with at least one choice of z re t- As mentioned in the 
last section, however, galaxy formation in low mass haloes can be 
highly stochastic. Here we use the observables predicted by our 
model to test whether or not if suffers from a TBTF problem. 

Fig.[5]shows the fraction of Milky Way analogues which host 
N or fewer satellites with Knax greater than a threshold value 
(given in the legend). In the top panel, as |Boylan-Kolchin et al.| 
(2dTT 20l2i did, we select the 12 subhaloes with the largest val¬ 
ues of Knax at 2 = 0. We immediately see the fraction of Milky 
Way analogues hosting no more than N = 3 satellites galaxies 
with Vmax > 30kms _1 is zero, in line with the Aquarius sim¬ 
ulations discussed above. However, assuming approximate y/N 
shot noise in the real Milky Way count, we can also note that 
the fraction of hosts having no more than N = 5 subhaloes with 
Knax > 30kms _1 is 23%. 

The fraction of ‘MW-compatible’ systems decreases at a given 
value of N if we lower the threshold Knax (red and cyan curves), 
and increases if we raise the threshold (blue curve). A recent study 
compared Local Group analogues in hydrodynamical and colli¬ 
sionless simulations from the same initial conditions jSawala et al.| 
|2014j > and found that baryonic processes associated with feedback 
systematically reduce Knax for a given subhalo by ~ 10% - 15%. 
As our simulation is collisionless, it is reasonable to take this ef¬ 
fect into account by assuming observed satellites with a measured 
Vmax = 30 km s -1 correspond to collisionless haloes with slightly 
higher Knax- This makes a significant difference - for example. 


stating the problem in terms of subhaloes with Knax > 35 kins -1 
yields a fraction of ‘MW-compatible’ hosts around 40%. 

The TBTF problem can be further alleviated by selecting sim¬ 
ulated satellites for comparison according to their observable prop¬ 
erties, rather than Knax- In the bottom panel of Fig. [5] we select 
the 12 satellites with the highest stellar mass and ask if this sample 
also suggests a TBTF problem. The figure shows that, even without 
considering the possibility that baryonic processes may reduce the 
values of K ma x, 10% of Milky Way analogues in our model have 
three or fewer satellites with K max < 30kins -1 , compatible with 
the real Milky Way. If systems with N = 5 or fewer satellites are 
considered comparable, then even a low Kn ax threshold 20 km s -1 
yields a compatible fraction of 10% (this is a very conservative 
threshold, since all 12 brightest satellites of the MW, except the 
LMC, SMC and Sagittarius, have Kirc < 20kms _1 at their half 
light radius). 

In summary, in addition to baryonic effects on satellite poten¬ 
tials and the uncertainty in the appropriate statistical error bar on 
the observed count of Milky Way satellites, the selection of simu¬ 
lated satellite galaxies by observables increases the fraction of sim¬ 
ulated ~ 10 12 Mg haloes that have abundances of ‘high Knax’ 
satellites similar to that of the Milky Way system. Therefore, al¬ 
though lowering the assumed host halo mass alleviates the TBTF 
problem jVera-Ciro et al.||2013| |Wang et al.||2012| |Cautun et al.| 
|2014| > this may not be necessary if a self-consistent galaxy forma¬ 
tion model, rather than a K ma x scaling relation, is used to identify 
the ‘brightest’ satellites. This approach may be crucial for formulat¬ 
ing realistic CDM-based theoretical predictions for MW-like satel¬ 
lite systems and their properties. 


3.3 “Dark” substructures 

The previous sections demonstrated that our model can reproduce 
the abundance, radial profile and Knax distribution of the satellites 
in the Milky Way, even though the abundance of subhaloes with 
high Knax is larger than the count of bright satellites observed with 
the same Knax- In this section, we will show more quantitatively 
the diverse properties of satellite galaxies at a given Knax- 

3.3.1 Knax vs. luminosity 

The top panel of Fig. [6] shows the relation between maximum cir¬ 
cular velocity and V-band magnitude for satellites in our 13 Milky 
Way analogues. The predicted median Knax is an increasing func¬ 
tion of V-band luminosity, but the scatter between Kn ax and magni¬ 
tude, My, is very large. Even haloes with Knax > 30 or 20 km s _1 
can host galaxies fainter than My = —9 (the approximate limits of 
completeness for current all-sky surveys of the Milky Way) or —5, 
suggesting that a significant number of 20 — 30kms _1 satellites 
are still to be discovered. Likewise, at a given magnitude, K max 
covers a wide range. For example, at My ~ —12, Knax varies 
from 40 km s _1 to < lOkms -1 . 

This scatter between K ma x and luminosity stems from the fact 
mentioned above that the star formation rate history of a dwarf 
galaxy depends not only on its final halo mass or Knax, but also 
on its dark matter mass assembly history |Li et al.|p010) . The bot¬ 
tom panel shows the relation between final stellar mass and Kn ax 
measured just before the time of infall. This relation is much tighter 
than the stellar mass vs. present-day Knax relation shown in the 
top panel. This is because the subhaloes hosting satellite galaxies 
are subject to tidal forces capable of stripping their mass. Galaxies 
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Figure 6. Top panel: V ma x vs. V-band luminosity for satellites of the Milky 
Way. Grey diamonds are our model predictions. Blue diamonds highlight 
satellites of the Milky Way analogue with a Vmax distribution consistent 
with the observed data over the whole magnitude range (blue curve in Fig. 
4). Thick and thin purple curves give the median and 16 to 84 percent range, 
respectively. Bottom panel: V^n ax i n faii vs. V-band luminosity for satellites 
of the Milky Way. The colour coding and line coding are the same as in the 
upper panel. Red circles with error bars are the observed MW satellites. Er¬ 
rors in magnitude are taken from Wolf et al.(2010), while those in Vkiax are 
estimated by Boylan-Kolchin et al. (2012). Green symbols show results for 
satellite galaxies in M31 as given by Tollerud et al. (2014). We do not in¬ 
clude the SMC and LMC in this figure. V ma x is ^ 50 — 60 km s~ 1 for the 
SMC jStanimirovic et al.|2004[|Harris & Zaritsky|2006| and > 80 km s —1 
for the LMC (Olsen et al. 2011). For Sagittarius, the appropriate value of 
Vrnax is hard to determine, because the satellite appeal's to be strongly per¬ 
turbed by its interaction with the central potential of the Galaxy. 


embedded deep inside the potential of the substructures are more 
condensed and thus more resistant to such stripping, which can 
lower their central stellar velocity dispersion without necessarily 
disturbing their surface brightness profile ( [Penarrubia et al,|2008| 
Cooper et al.|2010 ». Our model implies that in 1/10 of MW ana¬ 


logues ~ 1 — 2 relatively bright (My < —13) satellite galaxies 
may have extremely low values of V meLX (blue diamonds). In detail, 
this number will depend sensitively on the depth at which stars are 
embedded in their host potentials and on the details of tidal strip¬ 
ping. 


Fig ,|6]also shows the inferred relation between Vn ax and lumi¬ 
nosity for Milky Way and M31 satellites. V max is estimated in the 
same way in both cases {Tollerud et al.|20l4) . Most of these data 
lie below the median prediction of the model (thick blue curve), but 
well within the 16-84 percentile range. Interestingly, however, the 
most massive satellites in both observational datasets lie furthest 
below the predicted median relation. 


Figure 7. Upper panel: number of subhaloes with galaxies fainter than 
My = —9 as a function of V max . The black curve shows the median 
value for all 13 Milky Way analogues and the error bars show the 16 to 84 
percentile range. Red curves show the one system that is consistent with the 
observed Milky Way satellite Umax distribution, as shown in Fig. [4] Bottom 
panel: fraction of the subhaloes hosting galaxies fainter than My = —9 as 
a function of Umax. The colour coding is the same as in the upper panel. 


3.3.2 Detection fractions 

In the last section we showed that the scatter between Vn ax and 
luminosity is large, such that it is possible for relatively massive 
subhaloes to host a galaxy with luminosity below the current all¬ 
sky completeness limit. In this section we make an estimate of the 
fraction of undetected subhaloes as a function of Umax. For this 
purpose, we consider My = —9 to be the current limit of com¬ 
pleteness and further assume that this is a hard cut. (In practice 
this limit is also a function of surface brightness, and survey depth 
varies across the sky; e.g. |Koposov et al.|2008[ ) 

The top panel of Fig. [7] shows the average number of sub¬ 
haloes hosting galaxies fainter than My = —9 as a function of 
Vmax. At 20kms _1 we find a median of 18 subhaloes below this 
completeness limit. This corresponds to 45% of all the dark matter 
substructures of the same Vmax, as shown in the bottom panel. This 
fraction increases rapidly with decreasing 14n ax , such that around 
90% of subhaloes (~ 250) with Vmax > 10kins -1 are below 
the limit. This conclusion holds for the single Milky Way analogue 
with a Vmax distribution consistent with Milky Way observations 
(red curve in Fig.[4]l. 


4 CONCLUSIONS 

We have applied a semi-analytic galaxy formation model to cos¬ 
mological N-body simulations with sufficiently high resolution to 
study the properties of Milky Way analogues and their satellite 
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galaxies. Our two simulations give converged predictions for the 
dark halo mass of the Milky Way, despite a large difference in mass 
resolution (a factor of 60) and spatial resolution (a factor of 4). 
Given its stellar mass, the most probable value for the halo mass of 
the Milky Way according to our model is 1.4 x 10 12 Mq, with a 
la dispersion in the range 8.6 x 10 11 Mq to 3.1 x 10 12 Mq. This 
median value is consistent with measurements of the Milky Way 
halo mass using different dynamical tracers (see Wang et al.|2015[ 
for a summary). 

For any given halo mass, the probability of hosting a Milky 
Way analogue is rather low, with a maximum probability of ~ 20% 
in haloes of mass ~ 10 12 Mq. This value decreases rapidly with 
lower and higher host halo mass. This implies that if one wishes 
to simulate the formation of a Milky Way like system, halo mass 
selection alone is not sufficient: other factors such as environment 
and assembly history are important. We intend to explore the influ¬ 
ence of such factors on the formation of Milky Way analogues in 
future work. 

We used the COCO simulation to study the properties of satel¬ 
lites in Milky Way analogues. Our model is able to reproduce both 
the abundance of satellites as a function of V-band luminosity and 
the radial profile of the bright galaxies (My < —9). We have com¬ 
pared the Vinax distribution of the simulated satellites in COCO 
to inferences of K max from observations. We find one out of 13 
K max distributions for our Milky Way analogues overlaps the ob¬ 
served distribution. Flowever, we caution that the ‘observed’ values 
of K ma x are themselves dependent on calibration against collision¬ 
less TV-body simulations. 

With this caveat, a conservative upper limit on Knax for most 
known satellite galaxies (excluding the LMC, SMC and Sagittar¬ 
ius) is ~ 30kms -1 . The severity of the “too big to fail” problem 
depends on a number of systematic uncertainties, such as the error 
assumed on the count of known satellites, the effect of baryonic 
processes on K max (Sawala et al. 2014 and references therein) and 
the strength of the correlation between K ma x and luminosity. Our 
model shows that this last factor - the selection of the compari¬ 
son sample according to luminosity in the presence of a large scat¬ 
ter between luminosity and Knax - is important for predicting the 
number of MW satellites with high Knax. 

The large scatter between luminosity (stellar mass) and Knax 
is a consequence of the complex formation histories of dwarf galax¬ 
ies and environmental effects (such as tidal and ram-pressure strip¬ 
ping) acting on dark matter subhaloes. Subhaloes with high Kn ax 
can host very faint galaxies, and bright galaxies can be hosted by 
haloes of low Knax- A significant fraction of subhaloes in our 
model host galaxies below the approximate all-sky completeness 
limit for present surveys of My > —9. Around half of our sub¬ 
haloes with Knax > 20km s -1 host galaxies fainter than this 
limit. The continued discovery of fainter Milky Way companions 
by deeper surveys will be important for constraining even the mas¬ 
sive end of the subhalo mass function and associated tests of the 
CDM model. 
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APPENDIX A: NUMERICAL CONSIDERATIONS 
CONCERNING THE MAXIMUM CIRCULAR VELOCITY 
OF SATELLITES 

When haloes fall into larger systems, they orbit around the poten¬ 
tial centres as substructures, losing some of their mass to tidal strip¬ 
ping. Some substructures survive for a long time, while others may 
be disrupted shortly after infall. The finite particle mass of N-body 
simulations sets a minimum mass below which subhaloes cannot 
be resolved (in SUBFIND, subhaloes must contain at least 20 parti¬ 
cles). Since most subhaloes lose mass though tidal stripping, many 
will eventually fall below this limit and be Tost’ from the simula¬ 
tion. The lower the resolution of the simulation, the more rapidly 
subhaloes will be lost artificially as a result of this numerical limit. 
The identification of the surviving substructures is thus affected by 
numerical resolution. 

The MSII and the COCO simulations differ by a factor of 
60 in particle mass and by a factor of 4 in softening length. In 
Fig. |Al| we show the difference in the cumulative number of sub¬ 
structures as a function of maximum circular velocity in the two 
simulations. The difference is less than 20 percent at velocities, 
Knax > 25kms _1 , but increases very rapidly to 100 percent at 
15kins -1 . At 20kins -1 , MS-II resolves only half of the sub¬ 
haloes found in COCO. We note that MS-II adopts a higher value 
of os which results in more substructures. Hence, the curve in 
Fig. |A 1 | is a lower limit on the difference due to resolution. Studies 
of satellite kinematics with MSII should be restricted to the regime 
Knax > 30kms _1 , while most of the Milky Way satellites have 
Umax < 20 km s _1 , except for the LMC, SMC and Sagittarius (e.g 
|Boylan-Kolchin et al.|2012) . It is thus less reliable to use the MSII 
to study the properties of satellite of the Milky Way without a very 
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Figure Al. Upper panel: the ratio of the cumulative number of subhaloes as 
a function of Umax in the COCO and MS-II simulations. Bottom: the ratio 
of the cumulative number of satellites with and without resolved subhaloes 
as a function of Umax- 


careful treatment of these satellite galaxies whose subhaloes are 
lost due to resolution. 

The semianalytic galaxy formation model circumvents the res¬ 
olution limit by allowing satellite galaxies to survive as bound ob¬ 
jects for longer than their corresponding N-body subhaloes. We re¬ 
fer to these as “orphan” galaxies, because their parent subhalo has 
been lost. Their continued survival is determined by analytic cal¬ 
culations of inspiral towards the central galaxy under dynamical 
friction and the impact of tidal stripping on their stellar profile. The 
value of V'max for orphan galaxies is fixed to that of their parent 
subhalo at the time it was last resolved. However, since these ob¬ 
jects are suffering severe tidal stripping by definition, this estimate 
of Vmax is unlikely to be very accurate. 

The bottom panel of Fig. |A1| shows the fraction of satellite 
galaxies in our model whose subhaloes are resolved (i.e. are not or¬ 
phans) as a function Vjnax- We find more than 90 percent of satel¬ 
lite galaxies in COCO are in resolved subhaloes, and thus should 
have reliable Vjnax measurements. In the main body of the paper we 
therefore use results from the COCO simulation to study satellite 
galaxies. 








